Model Evaluation & ROC curve

################################################################################
# Model Evaluation & ROC curve                                                 #
# ref : machine learning with R, CH 10                                         #
# ref : https://www.r-bloggers.com/a-small-introduction-to-the-rocr-package/   #
# ref : https://www.r-bloggers.com/illustrated-guide-to-roc-and-auc/           #
################################################################################

################################################################################
# terms                                                                        #
# TP FN TF FP                                                                  #
# recall , precison , accuracy , error rate                                    #
# Sensitivity , Specifity                                                      #
# Kappa Stat                                                                   #
# F measure                                                                    #
# roc curve , auc                                                              #  
# cross validation                                                             #
################################################################################

################################################################################
# package & methods                                                            #
# CrossTable in gmodels package                                                # 
# confusionMatrix in caret package                                             #
# createFolds in caret package                                                 #
# ROCR package                                                                 #
################################################################################

#0. 기본적인 model evaluation 에 쓰이는 지표들은 위와 같다.
# 서비스에서 탐지하려는 것이 어떤 것인가에 따라, 모델의 성능을 평가하는 척도도 달라지기 마련인데, 
# 기본적으로 confusion matrix 의 네가지 값의 조합으로 만들어진다.

#1. prediction 결과 Data Frame 
# success : 실제 성공 여부 값 
# prob    : 성공할확률 예측값
# predict : threshold 0.6 을 기준으로, threshold 값보다 크면 성공이라고 판단한 최종 예측값 
head(result)

#1.1 기초적인 Accuracy , Precision, Recall 을 구해보자 ( + Sensitivity , Specificity )
evaluate_model <- function(df) {
  a <- filter(df, predict == 1 & success == 1)
  b <- filter(df, predict == 0 & success == 1)
  c <- filter(df, predict == 0 & success == 0)
  d <- filter(df, predict == 1 & success == 0)
  
  success_precision <- nrow(a) / (nrow(a) + nrow(d)) 
  success_recall    <- nrow(a) / (nrow(a) + nrow(b)) 
  fail_precision    <- nrow(c) / (nrow(b) + nrow(c)) 
  fail_recall       <- nrow(c) / (nrow(c) + nrow(d))  
  accuracy          <- (nrow(a) + nrow(c)) / nrow(df)
  
  evaluation_result <- as.data.frame(list(
    `success_precision` = success_precision, 
    `success_recall` = success_recall, 
    `fail_precision` = fail_precision,  
    `fail_recall` = fail_recall, 
    `accuracy` = accuracy
  ))
  
  return(evaluation_result)
}
print(evaluate_model(result))

# precision 과 recall 앞에 success / fail 이라는 prefix 가 있다. 
# 예측모델에서 주안점을 둬서, 탐지해 내야할 대상이 success 이냐 fail 이냐에 따라서 달라질 수 있다는 이야기인데, 
# 위의 confusion matrix 의 기준이 Actually Spam 이냐 Acutally Ham 이냐의 차이이다. 

#1.2  직접 구현을 해서 기초 평가 지표를 뽑았으나, 라이브러리를 사용해도 된다.

library(gmodels)
CrossTable(result$success, result$predict)

 
   Cell Contents
|-------------------------|
|                       N |
| Chi-square contribution |
|           N / Row Total |
|           N / Col Total |
|         N / Table Total |
|-------------------------|

 
Total Observations in Table:  88908 

 
               | result$predict 
result$success |         0 |         1 | Row Total | 
---------------|-----------|-----------|-----------|
             0 |     24344 |     13500 |     37844 | 
               |  3779.282 |  2908.251 |           | 
               |     0.643 |     0.357 |     0.426 | 
               |     0.630 |     0.269 |           | 
               |     0.274 |     0.152 |           | 
---------------|-----------|-----------|-----------|
             1 |     14320 |     36744 |     51064 | 
               |  2800.861 |  2155.332 |           | 
               |     0.280 |     0.720 |     0.574 | 
               |     0.370 |     0.731 |           | 
               |     0.161 |     0.413 |           | 
---------------|-----------|-----------|-----------|
  Column Total |     38664 |     50244 |     88908 | 
               |     0.435 |     0.565 |           | 
---------------|-----------|-----------|-----------|

 

library(caret)
confusionMatrix(result$predict, result$success)

Confusion Matrix and Statistics

          Reference
Prediction     0     1
         0 24344 14320
         1 13500 36744
                                         
               Accuracy : 0.6871         
                 95% CI : (0.684, 0.6901)
    No Information Rate : 0.5743         
    P-Value [Acc > NIR] : < 2.2e-16      
                                         
                  Kappa : 0.3618         
 Mcnemar's Test P-Value : 9.095e-07      
                                         
            Sensitivity : 0.6433         
            Specificity : 0.7196         
         Pos Pred Value : 0.6296         
         Neg Pred Value : 0.7313         
             Prevalence : 0.4257         
         Detection Rate : 0.2738         
   Detection Prevalence : 0.4349         
      Balanced Accuracy : 0.6814         
                                         
       'Positive' Class : 0              
                                         

# 1.3 개별적인 지표들은 이제 살펴봤다. 
# 그러니까, precision / recall / sensitivity / specificity 와 Precision 과 Recall 의 산술평균'스러운' Accuracy 까지
# 아직 안 본 것이 F-measure 인데, 위의 식을 보면 조화 평균(harmonic mean) 을 사용하고 있다. 

# F-measure 는 각각 표현되는 두개의 값 precision 과 recall 을 combine 하여 나타내기 위한 것이다.
# 산식은 간단하고, 높으면 대체적으로 precision 과 recall 이 모두 우수하다는 것을 말한다.

# 하지만 왜? 산술 평균 말고 조화 평균인가?
# ref : http://www.4four.us/article/2013/09/how-to-combine-two-values

# 위의 그림과 링크에서 설명되듯이, 
# precision , recall 중 한쪽이 지나치게 나쁜 경우 
# ex> 실제 데이터의 실패:성공 = 0.99:0.01 인 경우, 성공을 예측하는 경우 recall 은 굉장히 좋은 반면 precision 이 엉망일 수 있다.
# 이 경우, accuracy 를 통해 평가를 하면 도드라지지 않지만, f-measure 를 하면 값이 일정 수준 이상 올라갈 수가 없다.
# 주의 : 모든 경우에 accuracy / f-measure 를 들이대면 안된다.

# 2. ROC curve & auc 
# 모델의 민감도(Sensitivity)와 특이도(Specifity) 의 관계를 표현한 그래프.
# X 축 : Sensitivity == TPR 
# Y 축 : 1 - Speicifity == FPR 
# TPR 과 FPR 은 반비례 관계 
pred <- prediction(result$predict, result$success)
roc.perf = performance(pred, measure = "tpr", x.measure = "fpr")
plot(roc.perf, ylab = "sensitivity ( TPR )", xlab = "1 - specificity ( FPR )")

# auc = Area Under the ROC Curve 
# auc 는 ROC Curve 의 아래면적이 크면 좋다는 것을 의미한다.
library(ROCR)
auc.perf = performance(pred, measure = "auc")
library(pROC)
auc.pROC = auc(result$success, result$predict)
paste('ROCR :', auc.perf@y.values , ' pROC : ', auc.pROC)

[1] "ROCR : 0.681419991398182  pROC :  0.681419991398182"

# extra visualization 
plot_pred_type_distribution <- function(df, threshold) {
  v <- rep(NA, nrow(df))
  v <- ifelse(df$prob >= threshold & df$success == 1, "TN", v)
  v <- ifelse(df$prob >= threshold & df$success == 0, "FN", v)
  v <- ifelse(df$prob < threshold  & df$success == 1, "FP", v)
  v <- ifelse(df$prob < threshold  & df$success == 0, "TP", v)
  
  df$predict <- v
  
  print(table(df$predict))
  
  ggplot(data=df, aes(x=success, y=prob)) + 
    geom_violin(fill=rgb(1,1,1,alpha=0.6), color=NA) + 
    geom_jitter(aes(color=predict), alpha=0.6) +
    geom_hline(yintercept=threshold, color="red", alpha=0.6) +
    scale_color_discrete(name = "type") +
    labs(title=sprintf("Threshold at %.2f", threshold))
}
plot_pred_type_distribution(result, 0.6)

   FN    FP    TN    TP 
13500 14320 36744 24344 

# 만약 서비스에서 success case 보다 fail case 를 더 잘 찾기를 원한다고 할때, 
# 기본적인 ROC curve 는 이를 잘 판단해주지 못한다.
# 이런 경우, 각각의 가중치를 두어 수치화해볼 수 도 있다.
# 물론 둘다 잘하면 좋다. 하지만 그렇지 못하는 경우... 하나만 더 잘 찾아내기를 원한다면.. 가중치를 두어 수치화해볼 수 있다.
calculate_roc <- function(df, cost_of_fp, cost_of_fn, n=100) {
  tpr <- function(df, threshold) {
    sum(df$prob <= threshold & df$success == 0) / sum(df$success == 0)
  }
  
  fpr <- function(df, threshold) {
    sum(df$prob <= threshold & df$success == 1) / sum(df$success == 1)
  }
  
  cost <- function(df, threshold, cost_of_fp, cost_of_fn) {
    sum(df$prob >= threshold & df$success == 0) * cost_of_fn + 
      sum(df$prob < threshold & df$success == 1) * cost_of_fp
  }
  
  roc <- data.frame(threshold = seq(0,1,length.out=n), tpr=NA, fpr=NA)
  roc$tpr <- sapply(roc$threshold, function(th) tpr(df, th))
  roc$fpr <- sapply(roc$threshold, function(th) fpr(df, th))
  roc$cost <- sapply(roc$threshold, function(th) cost(df, th, cost_of_fp, cost_of_fn))
  
  return(roc)
}
roc <- calculate_roc(result, 3, 1)

plot_roc <- function(roc, threshold, cost_of_fp, cost_of_fn) {
  library(gridExtra)
  
  norm_vec <- function(v) (v - min(v))/diff(range(v))
  
  idx_threshold = which.min(abs(roc$threshold-threshold))
  
  col_ramp <- colorRampPalette(c("green","orange","red","black"))(100)
  col_by_cost <- col_ramp[ceiling(norm_vec(roc$cost)*99)+1]
  p_roc <- ggplot(roc, aes(fpr,tpr)) + 
    geom_line(color=rgb(0,0,1,alpha=0.3)) +
    geom_point(color=col_by_cost, size=4, alpha=0.5) +
    coord_fixed() +
    geom_line(aes(threshold,threshold), color=rgb(0,0,1,alpha=0.5)) +
    labs(title = sprintf("ROC")) + xlab("FPR") + ylab("TPR") +
    geom_hline(yintercept=roc[idx_threshold,"tpr"], alpha=0.5, linetype="dashed") +
    geom_vline(xintercept=roc[idx_threshold,"fpr"], alpha=0.5, linetype="dashed")
  
  p_cost <- ggplot(roc, aes(threshold, cost)) +
    geom_line(color=rgb(0,0,1,alpha=0.3)) +
    geom_point(color=col_by_cost, size=4, alpha=0.5) +
    labs(title = sprintf("cost function")) +
    geom_vline(xintercept=threshold, alpha=0.5, linetype="dashed")
  
  sub_title <- sprintf("threshold at %.2f - cost of FP = %d, cost of FN = %d", threshold, cost_of_fp, cost_of_fn)
  grid.arrange(p_roc, p_cost, ncol=2)
}
# threshold = 0.83 / FP 가중치 : 3 , FN 가중치 : 1 
plot_roc(roc, 0.83, 3, 1)

# 3. cross validation
# 아래의 그림처럼 특정 data set 을 n 개 의 bucket 으로 나누고, n 회 반복하여 train set 과 test set 을 변경하여 모델을 평가한다. 
# 이를 통해 cherry-picked 되는 케이스를 걸러낼 수 있다. 

# 3.1 data partition
success_dataset   <- filter(raw, success == 1)
train_set_success <- sample_frac(success_dataset, 0.2)
test_set_success  <- setdiff(success_dataset, train_set_success)
# train_set <- rbind(train_set_success, train_set_fail)
paste("total : ", nrow(success_dataset), " train : " , nrow(train_set_success), " test : ", nrow(test_set_success))

[1] "total :  80315  train :  16063  test :  64252"

# 3.2 bucket
library(caret)
folds <- createFolds(raw$success, k = 10)
str(folds)

List of 10
 $ Fold01: int [1:13019] 3 6 12 15 20 23 31 73 74 78 ...
 $ Fold02: int [1:13020] 9 18 22 47 59 68 76 111 113 154 ...
 $ Fold03: int [1:13020] 2 30 41 49 62 71 72 84 101 102 ...
 $ Fold04: int [1:13020] 28 34 36 37 45 46 48 51 52 58 ...
 $ Fold05: int [1:13020] 10 13 27 56 97 108 117 146 149 151 ...
 $ Fold06: int [1:13020] 17 24 40 50 54 55 61 63 65 70 ...
 $ Fold07: int [1:13020] 1 7 14 33 53 57 81 83 87 90 ...
 $ Fold08: int [1:13020] 11 19 42 43 44 66 69 82 106 107 ...
 $ Fold09: int [1:13019] 4 8 16 21 25 29 38 60 89 98 ...
 $ Fold10: int [1:13020] 5 26 32 35 39 64 77 88 91 96 ...

test <- raw[folds$Fold01, ]
table(raw$success)

    0     1 
49883 80315 

table(test$success)

   0    1 
4970 8049 

# 4. kappa stat
# Kappa Statistic compares the accuracy of the system to the accuracy of a random system
# 따라서, 이 값이 작을 수록 좋은 놈이다.
# 예를들어 정답인 놈들만을 가지고 kappa 를 구해보면, 아래와 같이 0이 된다. 
tp <- head(subset(result, success == 0 & predict == 0), 10000)
r <- rbind(tp, head(subset(result, success == 0 & predict == 1), 1))
r <- rbind(r, head(subset(result, success == 1 & predict == 0), 1))
r <- rbind(r, head(subset(result, success == 1 & predict == 1), 10000))
cm <- confusionMatrix(r$predict, r$success)
print(cm$overall)

      Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull AccuracyPValue  McnemarPValue 
     0.9999000      0.9998000      0.9996388      0.9999879      0.5000000      0.0000000      1.0000000 

# 그리고 틀린 놈들의 비중을 높여서 kappa 를 보면,
tp <- head(subset(result, success == 0 & predict == 0), 10)
r2 <- rbind(tp, head(subset(result, success == 0 & predict == 1), 1000))
r2 <- rbind(r2, head(subset(result, success == 1 & predict == 0), 1000))
r2 <- rbind(r2, head(subset(result, success == 1 & predict == 1), 10))
cm <- confusionMatrix(r2$predict, r$success)
print(cm$overall)

      Accuracy          Kappa  AccuracyLower  AccuracyUpper   AccuracyNull AccuracyPValue  McnemarPValue 
   0.009900990   -0.980198020    0.006057974    0.015250042    0.500000000    1.000000000    1.000000000 

# Kappa 통계량은 아래와 같은 기준을 가지고 생각하면 되고, 산식은 아래와 같다.
# Kappa Agreement
# < 0 Less than chance agreement
# 0.01–0.20 Slight agreement
# 0.21–0.40 Fair agreement
# 0.41–0.60 Moderate agreement
# 0.61–0.80 Substantial agreement
# 0.81–0.99 Almost perfect agreement

# 5. 적합한 measurement 는 어떤 것인가? 케이스 분석 

# 5.1 매우 조심스럽게 fail 케이스를 탐지하고, fail 이 예상되는 경우 보수적인 대책을 제시하는 것이 목표라면 
# precision of fail 이 제일 중요한 요소가 된다.  ( threshold = 0.1 선택 )
# 5.3 매우 공격적으로 fail 케이스를 탐지하고, FN 이 발생할 경우 risk 를 감당해낼 수 있다면,
# recall of fail 이 제일 중요한 요소가 된다. ( threshold > 0.7 선택 )
# 5.2 길 찾기 표지판에서 오른쪽 왼쪽을 맞추는 문제를 푸는 것이라면, 
# accuracy 혹은 f-measure 가 기준이 된다.

# 6. 정리 
# EDA 와 기본적인 cleaning 에 필요한 통계적 고찰과 
# 모델의 결과에 대한 measurement 방식과 선택에 대한 판단이 갖춰지면, 
# 적합한 기법을 찾아 적용하면 된다. 
# 모델을 튜닝하고 앙상블하는 과정도 중요하지만, 가장 중요한 것 세가지는 cleaning 과 feature eng. 그리고 적절한 판단법에 근거한 의사결정이다.